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Resume 

In the framework of noisy quantum homodyne tomography with efficiency parame- 
ter < ?7 < 1, we propose two estimators of a quantum state whose density matrix 
elements pm,n decrease Uke e-s(™+")'^^^ for fixed known B > Q and < r < 2. The 
first procedure estimates the matrix coefficients by a projection method on the pattern 
functions (that we introduce here for < ry < 1/2), the second procedure is a kernel es- 
timator of the associated Wigner function. We compute the convergence rates of these 
estimators, in L2 risk. 
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1 Introduction 



Experiments in quantum optics consist in creating, manipulating and measuring quantum 
states of light. The technique called quantum homodyne tomography allows to retrieve 
partial, noisy information from which the state is to be recovered : this is the subject of the 
present chapter. 

1.1 Quantum states 

Mathematically, the main concepts of quantum mechanics are formulated in the language of 
selfadjoint operators acting on Hilbert spaces. To every quantum system one can associate 
a complex Hilbert space TL whose vectors represent the wave functions of the system. These 
vectors are identified to projection operators, or pure states. In general, a state is a mixture 
of pure states described by a compact operator p onTi having the following properties : 

1. Selfadjoint : p = p* , where p* is the adjoint of p. 

2. Positive : p > 0, or equivalently {ip, pip) > for all ip £ TC. 

3. Trace one : tr(p) = 1. 

When Ti. is separable, endowed with a countable orthonormal basis, the operator p is iden- 
tified to a density matrix [pm,n]m,n£N- 

The positivity property implies that all the eigenvalues of p are nonegative and by the trace 
property, they sum up to one. In the case of the finite dimensional Hilbert space C^, the 
density matrix is simply a positive semi-definite d x d matrix of trace one. Our setup from 
now on will hen = L^(M), in which case we employ the orthonormal Fock basis made of 
the Hermite functions 



where H.m{x) := {—l)"^e^ m-th Hermite polynomial. Generalizations to 

higher dimensions are straightforward. 

To each state p corresponds a Wigner distribution Wp, which is defined via its Fourier 
transform in the way indicated by equation ^ : 



hm{x) := {1"'m\y/^)-'2H^{x)e-^ 



(1) 



Wp{u,v) : 




e-i^^i+MWp{q,p)dqdp := Tr(/9exp(-iuQ - ivP 



)) 



(2) 
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where Q and P are canonically conjugate observables (e.g. electric and magnetic fields) 
satisfying the commutation relation [Q, P] = i (we assume a choice of units such that 
^ = 1). It is easily checked that Wp is real- valued, has integral JJ^2 Wp{q,p)dqdp = 1 and 
uniform bound |Wp(q',p)| < ^. 

For any S M, the Wigner distribution allows one to easily recover the probability density 
X Ppix, (j)) of Q cos (j) + P sin (p by 

Pp{x,4>)=n[Wp]{x,<p), (3) 

where TZ is the Radon transform defined in equation 

/oo 
Wp{x cos (j) — t sin (j), X sin (j) + t cos (/))dt. (4) 

Moreover, the correspondence between p and Wp is one to one and isometric with respect 
to the L2 norms as in equation ([5]) : 

Il^pll2 := // \Wp{q,p)\'dqdp = ^\\p\\l := ^ (5) 

j,k=0 

From now on we denote by (•, •) and || • || the usual Euclidian scalar product and norm, while 
C(-) will denote positive constants depending on parameters given in the parentheses. 
We suppose that the unknown state belongs to the class TZ{B, r) ior B > and < r < 2 
defined by 

7l{B,r) := {p quantum state : \pm,n\ < exp(— i?(m + nY^"^)}. (6) 

For simplicity, we have chosen to express the results relative to a class which is the intersec- 
tion of the (positive) ball of radius 1 in some Banach space with the hyperplane tr{p) = 1. 
Another radius for the class would only change the constant C in front of the asymptotic 
rates of convergence that we will find. 

As it will be made precise in Propositions [T] and El quantum states in the class given in dHJ 
have fast decreasing and very smooth Wigner functions. From the physical point of view, 
the choice of such a class of Wigner functions seems to be quite reasonable considering that 
typical states p prepared in the laboratory do satisfy this type of condition. 
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1.2 Statistical model 

Let us describe the statistical model. Consider (Xi, . . . , $„) independent identi- 
cally distributed random variables with values in Rx [0, vr] and distribution Pp having density 
Pp{x,(f)) (given by JS]) with respect to :^A, A being the Lebesgue measure on M x [0, vr]. The 
aim is to recover the density matrix p and the Wigner function Wp from the observations. 
However, there is a slight complication. What we observe are not the variables (X^, but 
the noisy ones (Y^, <!>£), where 



(7) 



with ^£ a sequence of independent identically distributed standard Gaussians which are 
independent of all (Xj,<I>j). The detection efficiency parameter < r/ < 1 is known from 
the calibration of the apparatus and we denote by A^'' the centered Gaussian density of 
variance (1 — t/)/2, and its Fourier transform. Then the density pp of (Y^, <I>^) is given 
by the convolution of the density Pp{-/ with N'^ 



In the Fourier domain this relation becomes 



where !Fi denotes the Fourier transform with respect to the first variable. 

The theoretical foundation of quantum homodyne tomography was outlined in 



(8) 



291] and 



has inspired the first experiments determining the quantum state of a light field, initially 



with optical pulses in 



26 



25 



has been discussed or studied in 10^ ^, 201. ,^ 
methods have been studied in 



19l |. The reconstruction of the density from averages of data 
l|] for T] = 1 (no photon loss). Max-likelihood 
and procedure using adaptive tomographic 



kernels to minimize the variance has been proposed in Hj]. The estimation of the density 
matrix of a quantum state of light in case of efficiency parameter ^ < t] < I has been 
discussed in p, Q, Q] and considered in 23] via the pattern functions for the diagonal 



elements. 
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1.3 Outline of the results 

The goal of this chapter is to define estimators of both the density matrix and the Wigner 
function and to compare their performance in L2 risk. In order to compute estimation risks 
and to tune the underlying parameters, we define a realistic class of quantum states TZ{B, r), 
depending on parameters B > and < r < 2, in which the elements of the density matrix 
decrease rapidly. 

In Section [21 we prove that the fast decay of the elements of the density matrix implies both 
rapid decay of the Wigner function and of its Fourier transform, allowing us to translate 
the classes TZ{B, r) in terms of Wigner functions. 

In Section [Sj we give estimators of the density matrix p. The legend was somehow forged 
that no estimation of the matrix is possible when < r/ < 1/2. The physicists argue 
that their machines actually have high detection efficiency, around 0.8 ; it is nevertheless 
satisfying to be able to solve this problem in any noise condition. We give here the so- 
called pattern functions to use for estimating the density matrix in the noisy case with any 
value of T] between and 1. These pattern functions allow us to solve an inverse problem 
which becomes (severly) ill-posed when < < 1/2. In this case, we regularize the inverse 
problem and this introduces a smoothing parameter which we will choose in an optimal 
way. We compute the upper bounds for the rates achieved by our methods, with L2 risk 
measure. 

In Section HJ we study a kernel estimator of the Wigner function in L2 risk, over the 



same class of Wigner functions. It is a truncated version of the estimator in ^| and tuned 
accordingly. We compute upper bounds for the rates of convergence of this estimator in L2 
risk. 

To conclude, we may infer that the performances of both estimators are comparable. We 
obtain nearly polynomial rates for the case r = 2 and intermediate rates for < r < 2 (faster 
than any logarithm, but slower than any polynomial). It is convenient to have methods to 
estimate directly both representations of a quantum state. The estimator of the matrix p 
can be more easily projected on the space of proper quantum states. On the other hand, we 
may capture some features of the quantum states more easily on the Wigner function, for 
instance when this function has significant negative parts, the fact that the quantum state 
is non classical. 
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2 Decrease and smoothness of the Wigner distribution 

We recall that the Wigner distribution Wp was defined in the introduction. In the Fock 
basis, we can write Wp in terms of the density matrix [pm,n] as follows (see Leonhardt [l^ 
for the details). 

Wp{q,p)=Y^ 

Pm,n 

m,n 

where 

Wm^n{q,p) = ^ j e^'f^'kmiq - x)hn{q + x)dx. (9) 
It can be seen that Wm,n{(l-,p) = Wn^miQ, —p) and if m > n, 



i-irfniy (,.^p. 



vr V m! 



,1^ 

, m—n 



X 



V2iip-q)) Lr"(2g' + 2j>2) (10) 



thus, writing z := \fq^\-T? 



1 



/ n! \ 2 



Vn(^) := W^M^P)\ = ^— ( ) e-^'z'"-"|Lr"(2^')| (11 



vr ym! 

where L^{x) := (n!)^^e^x~"^^(e~^a;""'"'^) is the Laguerre polynomial of degree n and order 
a. Concerning the Fourier transforms, we also recall that 

wZr{q,p) = ^~'^2^" ^rn,n{l,l). (12) 

In this section we show how a decrease condition on the coefficients of the density matrix 
translates on the corresponding Wigner distribution. First the case r < 2 : 

Proposition 1. Assume that < r < 2 and that there exists B > such that, for all 

m > n, 

\Hm,n\ _ c 

Then for all (3 < B, there exists zq (depending explicitly on r,B,P, see proof) such that 



z := \J q^ ^p^ > zq implies 

\Wp{q,p)\<A{z)e-''''' (13) 

as well as 

Wpiq,p) < A(z/2)e-^(^/2)'- (^4^ 



where A[z) := ^ ( E e-^("^+")''^' + ^z^-^. 
If r = 2, the result is a little different : 

Proposition 2. Suppose that there exists B > such that, for all m > n, 

\n I < 
\ym,n\ ^ ^ 



Then there exists zq such that z := y'q^+p^ > zq implies 



\Wpiq,p)\ < A{z)e (15) 

as well as 

r(^/2) 



Wp{q,p) <A{z/2)e (i+Vb)'2^ ' ' (ig) 



for A{z) = i ( E e-s("^+") + — ^ 



Note that ^^^^^a < min(i?,l). Even when S is very large, we cannot hope to obtain a 

2 

faster decrease because is the decrease rate of the basis functions themselves (Lemma 

The proof of these propositions is defered to Appendix[5l More general results and converses 
are studied in Q|. Let us now state a few general utility lemmata. 

Lemma 1. Let y and w be two functions : [xo,+oo) (0,+C!o) such that y'{x) 0, 
w is hounded, satisfying the differential equations 

y"{x) = (j){x)y{x) 
w" {x) = il){x)w{x), 

with continuous (j){x) < ip{x), and initial conditions y{xo) = w{xo). Then for all x > xq, 
w{x) < y{x). 

Demonstration. Suppose that there exists xi > xq where w{xi) > y{xi). Then for some 
X2 G [3^0) we have w'{x2) > y'{x2) and w{x2) > y{x2). Consequently, for all x > X2, 
w"{x) — y"{x) > 0, and w'{x) — y'{x) > w'{x2) — y'{x2). When re — > cxd, liminf w'(x) > 
w'{x2) — y'{x2) > 0, which contradicts the boundedness of w. □ 



This lemma is used to prove a bound on the Laguerre functions. 
Lemma 2. For all m,n G N and s := \/m + n + 1, for all z > 0, 

ifi a o< z< s 

lm,n{^)<-{ . - - (17) 

^ [ e-(^-") if z>s. 

Demonstration. When z < s, the result follows from the uniform bound on Wigner functions 
obtained by applying the Cauchy-Schwarz inequality to ([9]). 

When z > s, L^(2z^) doesn't vanish and keeps the same sign as L"(2s^). Now, as it can 



be seen from 



271 . 5.1.2], the function w{z) := y/zlm,n{z) satisfies the differential equation 



w" = (4(2;^ — s^) + ^^r^)2:. On the other hand, y[z) := ^/slm,n{s)e '^^^ satisfies y" 



(4(z - sf - 2)y. When z>s, 



4(^ -sf-2< 4(^2 - s^) + (1^ 



from which we conclude with Lemma [T] that w{z) < y{z). □ 
Finally, a lemma to bound the tail of a series. 

Lemma 3. If v > and C > 0, there exists a zq such that z > zq implies 

m+n>z 

Demonstration. First notice that 

m-t-n-^z i^z 

When t > z and z is large enough, we have 

Jz Jz 

Cv 

which is what we needed to prove. □ 
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3 Density matrix estimation 



The aim of this part is to estimate the density matrix p in the Fock basis directely from the 
data (Yi, <I'j)i=i,...,n- We show that for < ?y < 1/2 it is still possible to estimate the density 
matrix with an error of estimation tending to as n tends to infinity (Theorem . In both 
cases (?7 > ^ and r/ < we construct an estimator of the density matrix {pj,k)j,k<N-i 
from a sample of QHT data. We give theoretical results for our estimator when the quantum 
state p is in the class of density matrix with decreasing elements defined in ([6]). 



3.1 Pattern functions 

The matrix elements pj^k of the state p in the Fock basis ^ can be expressed as kernel 
integrals : for all j, /c G N, 

1 



Pj,k = - I 1^ Pp{x,4>)fj,k{x)e 



-i{k-j 



(20) 



where fj^k = fk,j are bounded real functions called pattern functions in quantum homodyne 
literature. A concrete expression for their Fourier transform using Laguerre polynomials 
was found in [2^ : for j > A;, 



fkjit) = 27r2|t|VF,-fc(t,0) 



(21) 



where fkj denotes the Fourier transform of the Pattern function fkj- 

Let us state the lemmata which are used to prove upper bounds in Propositions [3], H] and El 

Lemma 4. There exist constants C2, Coo such that 

Ef+k=o WfkAl < C2Nf and Zf+k=o Wfkjlt < ^ooNf. 
This is a slight improvement over jl], Lemma 1]. 

Demonstration. By symmetry we can restrict the sum to j > k. For fixed k and j we have 



fk,j 


' = / 




^dt+ 1 


IkAt) 


2 

dt 




2 J\t\<2s 


J\t\>2s 
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(with s = y/k + j + 1). Because of Lemma [2j it is clear that the second integral is negligible 

2 

in front of the first one, which we simply bound by 4s fkj 

I I °° 
In view of (f^ . the main result in [18|| can be rewritten as follows : if /c > 35 and j — k > 24, 

then 



fk,j 



< 28887r2(i + l)2k~6. 



In consequence, for these values of k and j, 

2 



fk,j 



< C{jk~6 +j2k3). 



(22) 



(23) 



On the other hand, a classical bound on Laguerre polynomials found in [27|] yields that, for 

2 



fixed values of j — /c. 



< Ck3 , hence for all k > 35 and j — k < 24, 



< C{j2k3 + ke). 



(24) 



When A; < 35, we can use another result in 
of j — k, thus 



17| which gives 



fk, 



fk,j 



<Cj. 



< Ckej2 independently 



(25) 



Comparing (f23ll . (f2^ and ([25]l we see that when N is large enough, in the sum over < 
j, k < N, the terms k > 35, j — k > 24 dominate and (f23l) yields the first inequality. 
The second inequality is obtained by doing a similar computation, starting with H/j-./cHqq < 
fj^k and using ([22]l to bound 



j,k 



< C{j2k 6 + j2k6) 



when k>35 and j -k> 24. □ 
In the presence of noise, it is necessary to adapt the pattern functions as follows. From now 



on, we shall use the notation 



7 : = 



When i < < 1, we denote by f'^- the function 
which has the following Fourier transform : 

fljit) ■■= fkAtW ■ (26) 

When < r/ < i, we introduce a cut-off parameter 6 > and define /^'| via its Fourier 
transform : 

Then we compute bounds on these pattern functions. 
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Lemma 5. For 1 > rj > 1/2, there exist constants and C2d such that 



Ef+.=o a, , < C^2N^^e''^ and Ef+.=o fl 



Demonstration. The proof is similar to the previous one and we skip some details. Once 
again we assume j > k and write 



fV 



\t\<2s 



kjit) 



\t\>2s 



fkAt) 



(where s = \Jk + j + 1). Because of Lemma O the second integral is of the same order as 
the first one, which we bound by 



fk. 



e^^'"dt < C 



fk, 



2 , 



°° J\t\<2s 

In the sum we are considering the terms A; > 35 and j — k > 2A are dominant and, once 
again thanks to (f22]l . remembering that s = ^J j + k + 1, 



Jk,j 



hence the first inequality. 

The second inequality is, in the same fashion, based on 



fV 

Jk,j 



< 



fV 

Jk,j 



^ < c[j-ik~T2 I e^^^dt 

1 \ J\t\<2s J 



when > 35 and j — k > 24, and the bound on the sum readily follows. 



□ 



3.2 Estimation procedure 

For N := N{n) — > oo and 5 := S{n) — > 0, let us define our estimator of pj^k for < j + A; < 
A^- 1 by 

(28) 



where 



/J|;.(x)e-*(J-^')'^ if i < < 1 
(x)e-^(j-'=)<^ if < 77 < i. 
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using the pattern functions defined in ([26]l and (f27l) . We assume that the density matrix 
p belongs to the class TZ{B,r) defined in (l6|). In order to evaluate the performance of our 
estimators we take the L2 distance on the space of density matrices ||r — := tr(|r— pp) = 
X]jfc=o \ ~ PjA^ ■ consider the mean integrated square error (MISE) and split it into 
a troncature bias term &f(n), a regularization bias terms 6|(n) and a variance term (T^(n). 

N-l 2 

= \p3,k\^ + \E[plk\- Phk\ 

j+k>N j+k=0 
N-l 

+ E E\plk-E[p]j 

j+k=Q 

=: bl{n) + bl{n)+a\n). 

The following propositions give upper bounds for 6f(n), 62 (^) ^^'^ cr^(n) in the different 
cases rj = 1, 1/2 <r7<lorO<r/< 1/2 and r = 2or0<r<2. Their proofs are defered 
to Appendix [5l 

Proposition 3. Let pjj^ be the estimator defined by ( fl5|) . for < rj < 1, with (5 — > and 
N ^ 00 as n ^ CO, then for all B > and < r <2, 

sup < c^^2-r/2g-2B^'V^ (29) 

pen{B,r) 

where ci is a positive constant depending on B and r. 

Proposition 4. Let p'j^ be the estimator defined by / fl5|) . for < 77 < 1/2, with — > 00 
as n —> 00 and 1/5 > . In the case r = 2, for {3 := B /{I + VB)"^ there exists 02, while 
in the case < r < 2, for any {3 < B there exists ci and no such that for n > uq : 

sup 62(n)<C2iV2^^r-12g-J^-|(i-2V]v)^_ ^3^^ 
pe7^(B,r) 

Note that for 1/2 < < 1 we have 62 (^) = for all < r < 2 (/5^^ is unbiased). 



Proposition 5. For pj^ the estimator defined by ^28\) . 

sup cj^(n) < C3 if < < 1/2 (31) 



p&n{B,r) 



n 



sup a^{n) < 4 e^T^ if 1/2 < 77 < 1 (32) 



pe'R.{B,r) 



n 



^17 

sup o"^(n) < C3 if r/ = 1 (33) 
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where C3, are positive constants depending on rj. 

We measure the accuracy of pj^ by the maximal risk over the class TZ{B,r) 

limsup sup fn"^ E i y2 \p\- Pj,k \<Co. (34) 
n-00 pe7^{B,r) y.^^^gl ' J 

where Co is a positive constant and (/j^ is a sequence which tends to when n ^ 00 and it 
is the rate of convergence. Cases rj = 1 (no noise), ^ < rj < 1 (weak noise) and < ry < ^ 
(strong noise) are studied respectively in Theorems [U [2] and [3l 

Theorem 1. When rj = 1, the estimator defined in for the model where the 
unknown state belongs to the class TZ{B,r), satisfies the upper bound [34\ ) with 

'^n = log(n)3rn 

2 

obtained by taking N{n) := ^ ^ ' • 

Demonstration. With the proposed N{n) one checks that the bias (f29l) is smaller than the 
variance ([33]) which is bounded by a constant times log(n)3f n~^. □ 

Theorem 2. When ^ < rj < 1, the estimator defined in ^2S\) for the model where the 
unknown state belongs to the class TZ{B,r), satisfies the upper bound |^|) with 

_ r = 2, 

127+B B 

ip^ = log(n) 3(47+s) n 47+s 

with N(n) — ^°g(") ( 1 + 2 log(logn) \ 
wtin ly \n) — 2(47+5) ^ 3 log(n) J ■ 

-For0<r<2, 

= log(n)2-/2e-2^^W^^^ 

where N{n) is the solution of the equation 8-jN + 2BN^/'^ = log(n). 
In that case we have N{n) = ^ log(n) — log(n)''/^ + o(log(n)'"/^). 

Demonstration. When r = 2, the proposed N{n) ensures that the variance ((32]) is equivalent 

127+B B 

to the bias ([29]) . which is bounded by a constant times \og{n)'-^('^~<+B'^ n . 

When < r < 2, the proposed N{n) makes the variance ([32]) bounded by a constant times 

g-2B7V(n) / ^ which is smaller than the bias, the latter being bounded by a constant times 

jY(ra)2-r/2g-2BAr(n)'-/2^ 
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The asymptotic expansion of N{n) is a standard consequence of its definition by the equation 
87 + 2BiV'^/2 = log(n) . □ 

Theorem 3. When < i] < ^, the estimator defined in / fl5|) for the model where the 
unknown state belongs to the class TZ{B,r), satisfies the upper bound [34\ ) with 

= N2-r/2^-2BNrn 

where N and 5 are solutions of the system 



128)^ ~^ 2\5 

for arbitrary [3 < B in the case < r < 2 or 



^ + i(i-2./]V)2 + |j = log(n) 
^ + i(i - 2^/Ny - 2BN''/^ = (loglog(n))2 



(35) 



^ + i(i - 2VNf - f log(iV) = log(n) 



-25^ ^ 2^5 - ) -3 

IF ^ 2V<5 



— + 2\/iV)2 - 2BN - 31og(iV) = 



(36) 



with P := ^-^^^p in the case r = 2. 

Theses bounds are optimal in the sense that ([351) and ([361) are obtained by minimizing the 
sum of the bounds (l29]l. (1301) and (|3T1). 



Demonstration. We use the standard notations a{n) ~ b{n) if — > 1 and a(n) ^ b{n) if 
there exists a constant M < 00 such that it < < ^ for all n. 

M — b(n) — 

Let us first examine the case < r < 2. Remark that the left-hand term of the second 
equation in ([35]) is strictly negative when 1/5 = 2\/]V and increases to 00 with 1/5. This 
proves that the solution satisfies 1/5 > 2\/]V and that Proposition H] applies. Furthermore, 
if we suppose that is unbounded when n — > 00, then (up to taking a subsequence) by 
the first equation -^-jr^ ~ log(n) whereas, by subtracting the two, |} ~ log(n), which is 
contradictory. So 1/5 « \/iV and we deduce that log(n). Then ((30]) yields 

whereas ([3T]l gives 
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We see that the dominant term is the bound (f29ll on bf{n), hence the result. 
When r = 2, the same reasoning as above yields 1/6 > 2\fN , 1/6 ^ ^/N and N log(n). 
Then the right-hand side of and are of the same order as A^e~^^^, which is the 
bound JSni) on □ 



4 Wigner function estimation 



4.1 Kernel estimator 

We describe now the direct estimation method for the Wigner function. For the problem of 



estimating a probability density / : ^ M directly from data {X^, <I>^) wit^ 
we refer to the literature on X-ray tomography and PET, studied by [2^, 



1 density TZ[f] 



m, m, m and 



many other references therein. In the context of tomography of bounded objects with noisy 



observations 



13(1 solved the problem of estimating the borders of the object (the support) 



The estimation of a quadratic functional of the Wigner function has been treated in 



22|. 



For the problem of Wigner function estimation when no noise is present, we mention the 

n 

work by They use a kernel estimator and compute sharp minimax results over a class 
of Wigner functions characterised by their smoothness. In a more recent paper [4], Butucea, 
Gu^a and Artiles treated the noisy problem for the pointwise estimation of Wp ; however the 
functions needed to prove minimax optimality there do not belong to the class of Wigner 
functions that we consider here. 

In this chapter, as in {4], we modify the usual tomography kernel in order to take into 
account the additive noise on the observations and construct a kernel which performs 
both deconvolution and inverse Radon transform on our data, asymptotically. Let us define 
the estimator : 



w;i{q,p) 



1 



Trn 



K? ( q cos ^£ + p sin ^£ 



where < r/ < 1 is a fixed parameter, and the kernel is defined by 



K» = - 



\t\ 



exp(-mt)|t| ^^ ^ 1^ 

4vr m{t/^) ' ^ 2iV'?(t/^) 



(37) 



(38) 



and /i > tends to when n — > 00 in a proper way to be chosen later. For simplicity, let 
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us denote z = {q,p) and [z, (p] = (?cos(/> + psincj), then the estimator can be written : 



This is a one-step procedure for treating two successive inverse problems. The main diffe- 

n 

rence with the noiseless problem treated by is that the deconvolution is more 'difficult' 
than the inverse Radon transform. In the literature on inverse problems, this problem would 
be qualified as severely ill-posed, meaning that the noise is dramatically (exponentially) 
smooth and makes the estimation problem much harder. 



4.2 L2 risk estimation 

We establish next the rates of estimation of Wp from i.i.d. observations (1^, ^i), £ = 1, . . . ,n 
when the quality of estimation is measured in L2 distance. In the literature, L2 tomography 
is usually performed for boundedly supported functions, see [l^ and [2]]]. However, most 
Wigner function do not have a bounded support ! Instead, we use the fact that Wigner func- 
tions in the class Tl{B,r) decrease very fast and show that a properly truncated estimator 
attains the rates we may expect from the statistical problem of deconvolution in presence 
of tomography. Thus, we modify the estimator by truncating it over a disc with increasing 
radius, as n — > 00. Let us denote 

D{sn) ={z = {q,p) £ R2 : ||z|| < Sn} , 

where s„ — > 00 as n — > 00 will be defined in Theorem [H Let now 

W',;iz)=WlJz)Inis„){z). (39) 

From now on, we will denote for any function /, 

f\z)dz, 



2 



D{s„) 



and by D{sn) the complementary set of D{sn) in M^. Then, 



E 



E 



E 



2 

D{s„) 



+ WW, 



2 

D{Sn) 



+ 



E 



2 

D(s„) 



+ \\Wo\\4, 
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When replacing the L2 norm with the above restricted integral, the upper bound of the 
bias of the estimator is unchanged, whereas the variance part is infinitely larger than the 
deconvolution variance in As the bias is dominating over the variance in this setup, we 
can still choose a suitable sequence Sn so that the same bandwidth is optimal associated 
to the same optimal rate, provided that Wp decreases fast enough asymptotically. The 
following proposition gives upper bounds for the three components of the L2 risk uniformly 
over the class TZ{B,r). 

Proposition 6. Let (Y^, ^i), i = 1, . . . ,n be i.i.d. data coming from the model ^ and let 
Wj^ be an estimator (with h ^ as n ^ 00) of the underlying Wigner function Wp. We 
suppose Wp lies in the class Tl{B,r), with B > and < r < 2. Then, for Sn ^ 00 as 
n ^ 00 and n large enough, 



< Cis. 



sup 

peniB,r) 



E[W;i] - Wp 



sup E 

pen{B,r) 



wl 

h,n 



E 



wJ! 

h,n 



D{s„) 

2 

DiSn) 



< Ca^exp 
nh 
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/l2 



where (3 < B is defined in Proposition{^ for < r < 2 and j3 = B/{1 + vB) for r = 2, 
7 = (1 — ri)/{Arj) > 0, Ci, C2, C3 are positive constants, Ci, C2, depending on (3, B, r and 
C3 depending only on rj. 

We measure the accuracy of W^'* by the maximal risk over the class TZ{B, r) 



lim sup sup E 

n^oo peTl{B,r) 



Wr. 



(40) 



where C is a positive constant and 99^ is a sequence which tends to when n — > 00 and it 
is the rate of convergence. 

In the following Theorem we see the phenomenon which was noticed already : deconvolution 
with Gaussian type noise is a much harder problem than inverse Radon transform (the 
tomography part). 

Theorem 4. Let B>0, 0<r<2 and {Y^^^i), i = 1, . . . ,n be i.i.d. data coming from 
the model f^. Then W^'* defined in ^39^) with kernel in liS8\} satisfies the upper bound 
{"JW with 
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- Forr = 2, put (5 = B/{1 + VB^ 



167+3g 



ipf^ = (log n) 8^+2/5 n , 



with Sn = {h) ^ andh= [^j:^ logn + log(logn) 
For < r < 2 and (3 < B defined in Proposition CI 



-1/2 



exp 



where Sn = ^/h and h is the solution of the equation 

h -p- = log n - (log log n) . 



Sketch of proof of the upper bounds. By Proposition [6j we get 

2^ 



sup E 

Wpen{B,r) 



< Cisi°-^'^e-"'< + C2h''-'^exp 



2/3 



27^ 

/l2 / ■ 



nh 

A1 + A2 + A3 



For < r < 2 and by taking derivatives with respect to h and s„, we obtain that the 
optimal choice verifies the following equations : 

2/3< + P = log(n)+log(/i4(4-'^)) 
^ + 5 = log(n)+log(/i2^-^.;^). 

We notice therefore that A2 is dominating over A3, which is dominating over Ai. The 
proposed (s„, h) ensure that the term A2 is still the dominating term and gives the rate of 
convergence. 

The case r = 2 is treated similarly, by taking derivatives we notice that the term A2 and 
the term A3 are of the same order and that the term Ai is smaller than the others. □ 
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5 Appendix 

5.1 Proof of Proposition [l] 

Let 0(2:) := {z — v^z''/^)^ — 1. Since r < 2, for z larger than a certain zq (which depends 

/ \ 2/r 

only on /3, B and r), it is true that (p{z) > ( j z^. It follows that 



(41) 



If m + n < (/)(z), then s < y^l + </>(z) and z - s > z - y^l + 0(z) = V^z'^/^^ By ([T7l), this 
means that lm,n{z) < ^e'^^"^ . So 

X] |Pm,n|/m,n(2) < Ae"^^'' (42) 
r?i+n<(/i(z) 

for 4 •= i V e~-^(™+")''''^ 

On the other hand, using Lemma [3] with v := r/2, if (p{z) > zq, 

V |p™,nK„.,n(^) < ^</>(z)2-/2e-^^W^^^ 



m+n>4i{z) 



< ^z^-^e"^^'' (43) 



by (flTll and (l4T]l . Combining ([42]) and ([43]) yields the announced result. The bound on Wp 
is then a direct consequence of ([T2]) . 



5.2 Proof of Proposition [2] 

Let 4>{z) := ez'^ - 1, where 9 := is the solution in (0, 1) of (1 - VO)'^ = Bd. 

When m + n < (j){z), then s < \fQz and z- s>z(\-^) = ^fmz. By ([T7D, this means 
that lm,n{z) < ie-™^'. So 

\Pm,n\lm,n{z) < ^e"^^^' (44) 

m+n<<j){z) 

for 4 ■= i V p-B(m+n) 

On the other hand, by Lemma [U if 4>{z) > zq, 

2 



\Pm,nK,n{z) < — (/.(z)e-^<^(^) 



m+n><f>{z) 



20e 



B 



ttB 



2 
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by (fT7|l and (|4T]) . Combining ([441) and ([451) yields the announced result. The bound on 
is then a direct consequence of (fT2]l . 



5.3 Proof of Proposition [3] 

By ([6]) the term (n) can be bounded as follows 

blin) = E exp(-2i3(j+A;r/2). 

j+k>N j+k>N 

Compare to the double integral and change to polar coordinates to get 

blin) < ciN^-'/'^exp{-2BN'/^). 

5.4 Proof of Proposition [4] 

To study the term b2{n), we denote 

^i[Pp{-mit) ■■= Ep[e''''\^ = <P] = Wp{t cos cP,t sin 
the Fourier transform with respect to the first variable. 



Y 



Y 



J, 



IT 



'i{j-k)4, I fV,i 



fj,'k iy)Vvp'liyVv\<P)dyd(p 



Jo ^ 



IT 



^-i{j-k)(t> 



1 

1 



-^i,-k)4>_ / fflit)T,[^pl{-^\<Pmdtd<l^ 



2tt 



\t\<i/& 



fj,k{t)e""'Ti[pp{-m{t)N\t)dtd<j). 



As N'^{t) = e-^^ and by using the Cauchy-Schwarz inequality we have 



vr 7o 27r 



fj,k{t)J'i\pp{-\ct>)]{t)dtd<t> 

\ 2 



< 



TT 



2tt 



|t|>l/5 



|t|>l/5 
f j,k{t)Wp{t COS 4), tsm. 



dt 



If 1/(5 > 2^/N > 2s with s = j + A; + 1, then whenever i > 1/5 we get by Lemma [2] 

\Jj,k{t)\ = 7r2|t|/,.fc(t/2) 

< 7r|t|e"3(l*l-2^)'. 
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On the other hand, by Propositions [T] and [2] we have 

\Wp{t cos (j),t sin < A(^)e~'^(^)' 



for P :-- 



in the case r = 2, or for arbitrary (3 < B and t large enough in the case 



< r < 2. In both cases ^ is a polynom of degree 4 — r. We deduce the inequahty 



by Lemma 8 in [5||, hence 



< C(-)12-4rg-i(i-2^)2-/32i-'-{|)'- 

6 



b2{nf < CiV2(l)12-4rg-i(i-2v/]V)2-/32i-(ir 



which covers both cases in the proposition. 



5.5 Proof of Proposition [5] 

Let us write crl^{n) := E plj^ - 



. We bound it by 



E 



1 " / Y 



n 



n 

< 

n 



2 



(46) 



Proof of (1311) For < r/ < 1/2, let us denote by Ks the function with the following 
Fourier transform Ks{t) = l{\t\ < \)e<^ , then /J^f = fjAt)Ks{t) and we have 



alj,{n) < ^E 



1 

n 



< -E 

n 



< -E 

n 



fj,k * Ks{ 



Y 



I fj,k{ms{^-t)dt 
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By using the Cauchy-Schwarz inequality 

1 



< 



\hkiWdtE 



dt 



1 



\fj^,{t)\'dtE— I Ks{u)e '"v^ 



1 2 /-i/^ 



2tt 



du 



Then, 



rnr ^ ^ \ — ri 

j+k=o ' 

By Lemma a we have J2f+k=o WfjMll < C2N^'^/^ thus 



n7r(l — rj) 

Proof of ([32]) and ([33]) By ([28]), for 1/2 < r? < 1, 



< 
< 



n 
1 

nvr 
1 

rnr 



-j(j-fc)$ 




/,'!fc(y) V^PpiVvy\4>)dyd(j) 



For 1/2 < ?7 < 1, by Lemma [5] 



nvr 



For 7/ = 1, by Lemma [6] 



n 




\fj,k{x)\^Pp{x,(t})dxd(t) 



hence by Lemma H] 



< — ll/j.fcllo 



a (n) < C . 

n 
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5.6 Proof of Proposition [6] 

It is easy to see that 



mi\\ H = Wp{w)i{\\w\\ < i/h). 

We have, for n large enough s„ > zq and by (fTSll 

||T^p|||j(,„) < C(i?,r) /" \\z\f-^^ eM-m41dz 

J ||2:||>5n 

/'27r /"CXD 

< C{B,r) / / t^-'^'' exp{-2(3f)dtdcj) 

JO J s„ 



where (3 < B and for n large enough in the case < r < 2, respectively (3 = B/{\ + 
in the case r = 2. Now we write for the L2 bias of our estimator : 



\E[WJI] -W-^ 



w, 



P\\2 



1 



< 



(2^)2 

C'^{B,r 



Wp{w) I{\\w\\ > l/h)dw 



\w 



|2(4-r) -2I 



dw 



w\\>l/h 



e 



by the assumption on our class and ([T4|, for < r < 2. The case r = 2 is similar. 
As for the variance of our estimator : 



V 



E 



2 

D{Sn) 



TT^n 



E 



E 



Y 



(47) 



On the one hand, by using two-dimensional Plancherel formula and the Fourier transform 
shown above, we get : 

2 



E 



Y 



\Wp{w)\^dw < T?. 



(48) 
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In the last inequality we have used the fact that ||Wp||2 = Tr(p^) < 1 where p is the density 
matrix corresponding to the Wigner function Wp. On the other hand, the dominant term 
in the variance will be given by 

2 



E 



Y 



D{s„) 



[K'^iiz, <P] - y/^)Y dzp^iy, ^)dydcP 

i^hi^))'^ VvP^iii^^ 'PI - u)VVi (t))dudzd(t) 
J 

[Kl{u)f [ r pp{-,<l)) * NN^iiz,^] - u)d(t>dzdu 

J D{Sn) Jo 



< M{r,)^sl j {Kl{u)fdu, 
using Lemma El below and the constant M(ry) > depending only on rj, defined therein. 



Indeed, let us note that y^Pp(-y^, </>) is the density of = X + — r])/{2r])e and 

let us call A^A^'' the Gaussian density of the noise as normalized in this last equation. 
Let us first compute, by Plancherel formula, H-K'^Hl and get 



dt 



1. f'^'' 



t^exp ( t^^^^ 1 dt 



exp (1;^) (1 + o{l)), ash^O. 



Anh 1 — r] 

We replace in the second order moment, then as /i ^ 



E 



Y 



D(s„). 



< 
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(49) 



The result about the variance of the estimator is obtained from (l47l)- (149)1. 



Lemma 6. For every p £ 1Z{B, r) and < < 1, we have that the corresponding probability 
density pp satisfies 

f'TT 

< / Ppi;C^)*NN'^{x)d(l)<M{7]), 

Jo 

< / Pp{x,<p)d(l)<C 
Jo 
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for all X G M eventually depending on (p, where M(r]) > is a constant depending only on 
fixed rj and C > 0. 



Demonstration. Indeed, using inverse Fourier transform and the fact that 
get : 

/"TT 



< 



< c{ri) j j \Wp{t cos (j),t 



sm ( 



exp 



4r] 



dtd(f) 



< c{r}) 



w 



Wp{w] 



exp 



4r] 



dw < M{r)), 



where c{r]), M(r]) are positive constants depending only on 77 G (0, 1). 



< 1 we 



□ 
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